The many-body exchange-correlation hole at metal surfaces 
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We present a detailed study of the coupling-constant-averaged exchange-correlation hole density 
at a jellium surface, which we obtain in the random-phase approximation (RPA) of many-body 
theory. We report contour plots of the exchange-only and exchange-correlation hole densities, the 
integration of the exchange-correlation hole density over the surface plane, the on-top correlation 
hole, and the energy density. We find that the on-top correlation hole is accurately described by local 
and semilocal density-functional approximations. We also find that for electrons that are localized 
far outside the surface the main part of the corresponding exchange-correlation hole is localized at 
the image plane. 

PACS numbers: 71.10.Ca,71.15.Mb,71.45.Gm 



I. INTRODUCTION 

The exchange-correlation (xc) energy of a many- 
electron system is the only density functional that has to 
be approximated in the Kohn-Sham (KS) formalism of 
density- functional theory (DFT)ii It is formally defined 
by the following equation derived from the Hellmann- 
Feynman theorem: 2 
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where n(r) is the density of a spin-unpolarized system 
of N electrons, U[n] = (1/2) / drn(r)n(r')/\r - r'| is the 
Hartree energy, and p 2 (r' , r) is the reduced two-particle 
density matrix 
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Here, 9 {r±cr±, tnctn) is the antisymmetric wavefunc- 
tion that yields the density n(r) and minimizes the ex- 
pectation value of f + \V ee , where f = - Y^Li Vf/2 
and V ee = \ J2i J2j^i \ri- rj \ are tne ki net i c energy and 
the electron-electron interaction operators. Eq. ([2]) shows 
that p 2 (r', r)dr'dr is the joint probability of finding an 
electron of arbitrary spin in dr' at r' and an electron of 
arbitrary spin in dr at r, assuming that the Coulomb in- 
teraction is A/|r — r'|. In the case of noninteracting KS 
electrons (i.e., A = 0), P2 = °( r '> r ) i s the exchange- only 
reduced two-particle density matrix that is expressible in 
terms of KS orbitals. (Unless otherwise stated, atomic 
units are used throughout, i.e., e 2 = h = m e = 1.) 

Hence, the xc energy can be expressed as the elec- 
trostatic interaction between individual electrons and 



the corresponding (and sorrounding) coupling-constant- 
averaged xc hole density fi xc ([n]; r,r'), as follows 

ip ri /"j / \ 1 f , f , ,n(r)n xc ([n];r,r') 
E xc [n] = J dr e xc (r) = - J dr J dr' , 
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where [see Eqs. (P and ©]: 
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and e xc (r) is the xc energy density. The xc hole density 
n IC ([n]; r, r') is the result of three effects: self-interaction 
correction to the Hartree approximation, Pauli exclusion 
principle, and the electron correlation due to Coulomb 
repulsion between electrons. 

The adiabatic-connection fluctuation-dissipation the- 
orem provides an elegant path to the exact coupling- 
constant-averaged xc hole densityj 3 -^^ which can be 
written as follows^ 
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where x A ( r : r '; u ) is the density-response function of the 
interacting system at coupling strength A and satisfies, in 
the framework of time-dependent density-functional the- 
ory (TDDFT), the following exact Dyson- type equation 8 

XA(r,r';w) = xo(r,r';cj) + J dr x dr 2 Xo(r,ri;u;) 



+ f xc ,x[n}(ri,r 2 ;ij) > x\(r 2 ,r';w). (6) 
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Here, xo(r, r'\u) is the density-response function 
of non-interacting KS electrons (which is exactly 
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known in terms of KS orbitals 9 ) and f xc ,\[n](r,r';uj) 
is the Fourier transform with respect to time 
[fxc,x[n](r,r';Lo) = dte*"/ BC , A [n](r, t, r>, 0)] of the 
unknown A-dependent xc kernel, formally defined by 



f xcA [n](r,t,r',t')= Sr - 
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where v xc [n](r,t) is the exact time-dependent xc poten- 
tial of TDDFT. When f xcA [n](r,r';uj) is taken to be 
zero, Eq. |6]) reduces to the random phase approxima- 
tion (RPA). If the interacting density response function 
X\( r i r 'i 0J ) is replaced by the noninteracting KS density- 
response function x ( r : r '; w ), then Eq. {5) yields the 
exchange-only hole density. 

The scaling relation of the correlation hole density at 
coupling constant A 10 ' 11 leads to the following equation 
for the coupling-constant-averaged correlation hole den- 
sity: 

n c ([n];r,r') = / dX (-) 3 t#(K/a], -r, -r'), (8) 
J Q w ww 

where < jo << 1 is a fixed constant, and n 7 (r) = 
7 3 n(7r) is a uniformly-scaled density— Eq. ([8} shows 
that the whole many-body problem is equivalent to the 
knowledge of the universal correlation hole density at a 
small, fixed coupling strength w. 

There is a "Jacob's ladder"- 3 classification (in RPA 
and beyond RPA) of nonempirical approximations to the 
angle-averaged xc hole density 

n xc ([n];r,u) = ^- J dfl n xc ([n}; r, r'), (9) 

where d£l is the differential solid angle around the di- 
rection of u = r' — r. The simplest rung of the lad- 
der is the local spin density approximation (LSDA) of 
the xc hole density n xc (n^,n^u) that has as ingredi- 
ents only the spin densities. (For the RPA-based LSDA 
xc hole and for the LSDA xc hole, see Refs. EH1 
and Refs. Il3fl6lfl7l . respectively.) The next rung is 
the generalized gradient approximation (GGA) xc hole 
density h xc (n^, nj,, Vn-j-, Vn.^, u). (See Ref. 14 for the 
smoothed GGA exchange hole model, Ref. [l8| for the 
PBE-GGAI2 correlation hole, and Ref. QJ for the RPA- 
based GGA hole model. For a GGA xc hole constructed 
for solids, see Ref. |20|.) The third rung on this lad- 
der is the non-empirical meta-GGA xc hole density 21 
n xc (n-\ , nj,, Vnj, Vn^, Tf,Tj_, u) that depends on spin den- 
sities and their gradients, as well as the positive KS ki- 
netic energy densities tj and tj , and that was constructed 
to satisfy many exact constraints. (For an RPA-based 
meta-GGA xc hole model, see also Ref. l2lh 

Jellium is a simple model of a simple metal, in which 
the ion cores are replaced by a uniform positive back- 
ground of density n = 3/47rr 3 = kp/Zir 2 and the va- 
lence electrons in the spin-unpolarized bulk neutralize 
this background. r s is the bulk density parameter and 



kp is the magnitude of the bulk Fermi wavevector. At 
a jellium surface, the plane z = separates the uniform 
positive background (z > 0) from the vacuum (z < 0), 
and the electrons can leak out into the vacuum. This 
electron system is translationally invariant in the plane 
of the surface. 

The exchange hole at a jellium surface was studied in 
Ref. |22 | ( using a finite linear-potential model^ 3 .) , and in 
Refs. 124051 (using the infinite barrier model (IBM)2&). 
The behavior of the xc hole at a jellium surface was in- 
vestigated at the RPA level using IBM orbitals. 27 Hence, 
existing calculations of the exchange-only and xc hole 
at a jellium surface invoke either a finite linear-potential 
model or the IBM for the description of single-particle 
orbitals. An exception is a self-consistent calculation of 
the RPA xc hole density reported briefly in Refs. [28| and 
[29l in which accurate LSDA single-particle orbitals were 
employed. 

In this paper, we present extensive self-consistent cal- 
culations of the exact-exchange hole and the RPA xc hole 
at a jellium surface. We report contour plots of the cor- 
responding hole densities, the integration of the xc hole 
density over the surface plane, and the on-top correla- 
tion hole. We find that the on-top RPA correlation hole 
n c ([n];r, r) is accurately described by the on-top RPA- 
based LSDA hole, in accord with the work of Perdew et 

n i 5.30,31 



II. THE EXACT-EXCHANGE HOLE AND THE 
RPA XC HOLE AT A JELLIUM SURFACE 

Let us consider a jellium surface with the surface plane 
at z — 0. Using its translational invariance in a plane per- 
pendicular to the z axis, the coupling-constant-averaged 
xc hole density of Eq. ([5]) can be written as follows 2 ^ 
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where r = |r|| — rji|, and q|| is a two-dimensional (2D) 

wavevector. X X (l\\ 1 z i z ' 1 LJ ) represents the 2D Fourier 
transform of the interacting density response function of 
Eq. ([6]), which in the RPA is obtained by neglecting the xc 
kernel f xc . The exact-exchange hole density is obtained 
by simply replacing in Eq. (fTU|) x A (<Z|| , z, z', oj) by the cor- 
responding KS noninteracting density response function 
X°(Q\\>z,z',uj). 

For the evaluation of Eq. (fTTJ]) . we follow the method 
described in Ref.0- We consider a jellium slab, and we as- 
sume that the electron density n{z) vanishes at a distance 
zo = 2\p (Xp — 2ir/kp is the bulk Fermi wavelength) 
from either jellium edge^ We expand the single-particle 
wave functions entering the evaluation of x (Q\\> z > z' t u>) 
in a sine Fourier representation, and the density-response 
functions x°(<Z|| ? z , z', oj) and x\{z, z';q\\,u>) in a double- 
cosine Fourier representation. We also expand the Dirac 
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delta function entering Eq. (fT0|) in a double-cosine rep- 
resentation (see Eq. (A2) of Ref. 0). We take all the 
occupied and unoccupied single-particle orbitals and en- 
ergies to be the LSDA eigenfunctions and eigenvalues of a 
KS Hamiltonian, as obtained by using the Perdew-Wang 
parametrizatior>2ii of the Ceperley- Alder xc energy of the 
uniform electron gasi^i ^ - 

In the calculations presented below, we have consid- q ^ 

ered jellium slabs with several bulk parameters r a and a 
thickness a = 2.23 Xp for the positive background. For _1 " 

r s = 2.07, such slab corresponds to about four atomic lay- _^ 5 _ 

ers of Al(100) and it was used in the wavevector analysis 
of the RPA2£ and beyond- RPA2£i2£ xc surface energy. ~ 2 2 

In Figs. [T]and [2J we show contour plots for the exact- 
exchange hole density and the self-consistent RPA xc hole 
density, respectively. In the bulk, both the exchange- 2 i— 

only hole and the xc hole are spherical and the xc hole 
is more localized, as in the case of a uniform electron 
gas. Near the surface, both the exchange-only hole and 1 - 

the xc hole happen to be distorted, the center of gravity 
being closer to the surface when correlation is included. 
For an electron that is localized far outside the surface, - 

the corresponding exchange-only hole and xc hole remain 1-1 ^ ^ 
localized near the surface; Figs. [1] and [2] show that the 
introduction of correlation results in a flatter hole, which _1 " 

in the case of an electron that is infinitely far from the 
surface becomes completely localized at a plane parallel 
to the surface. This is the image plane. We recall that ~ 2 2 15 1 Q5 Q Q5 1 15 

the RPA xc hole density is exact in the limit of large sep- z ^ 
arations (where u = |r — r'| — > 00), and yields therefore ' 
the exact location of the image plane. 2 

The integration of the xc hole density over the whole 
surface plane, 15 
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represents a quantity of interest for a variety of the- 
oretical and experimental situations (see for example -0.5 - 
Refs. HII). Below we show that b xc ([n]; z, z') repre- 
sents a suitable quantity to describe the behavior of the 
xc hole corresponding to a given electron located at an -1.5 - 

arbitrary distance from the surface. In Fig. [3J we plot _ 2 I 1 1 1 i 1 1 . 

this quantity, versus z' , for r s = 2.07 and a given elec- " 2 -1-5 -1 -0.5 0.5 1 1.5 

tron located at z = 0.5Af, z = 0, z = — 0.5Af, and zA, F 
z = — 1.5Ajr. We see from this figure that (i) correla- 
tion damps out the oscillations that the exchange hole 2 I ' ' ' ! ' ' ' 
exhibits in the bulk part of the surface, and (ii) in the 
case of a given electron located far from the surface into 
the vacuum the main part of the exchange-only and the 
xc hole is found to be near the surface (see also Figs. [1] 0.5 
and [2]) , although the exchange-only hole appears to be 
much more delocalized with a considerable weight within 
the bulk. -0.5 

Let us now focus on the on-top xc hole. The LSDA ac- 
curately accounts for short wavelength contributions to 
the xc energyj^S, thus, all the nonempirical approxima- -1.5 - 

tions of the xc hole have been constructed to recover the 2 | , , , j , , , 

LSDA on-top xc hole n^ DA (r,r). The slowly-varying -2 -1.5 -1 -0.5 0.5 1 1.5 
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FIG. 1: Contour plots of the exchange hole density 
n x {ru, z, z') for several fixed values of the electron posi- 
tion: 2 = 0.5Af (inside the bulk), z — (on the surface), 
7. — — f) fiAc fin thn variniml and 7. — —1 KAzr ('far mitsirle the 
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FIG. 2: Contour plots of the RPA coupling-constant- 
averaged xc hole density n m (rn, z, z') for several fixed values of 
the electron position: z = 0.5Af (inside the bulk), z — (on 
the snrfaceV 7. = — D.KAci (\r\ the variniml anH 7. = — I JSAc 



FIG. 3: b xc (z, z') of Eq. (JTTJ versus z' /\f for the same posi- 
tions of the electron as in Figs. [T]and[2] The bulk parameter 
is r s = 2.07 and the jellium surface is at z — 0. 
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FIG. 4: On-top coupling-constant-averaged correlation hole 
n c (r, r) at a jellium surface. Also shown is b c (z,z) of Eq. 

The bulk parameter is r s — 2.07 and the jellium surface 
is at z = 0. 



FIG. 5: The correlation hole b c (z,z) of an electron at po- 
sition z = — 0.5Af for several values of the bulk parameter 
r s = 2.07, 3, 4, 5, and 6. The jellium surface is at z = 0. 



electron gas was treated within RPA by Langreth and 
Perdew 5 -. For a spin-unpolarized system, the gradient 
correction to the LSDA on-top correlaton hole density 
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In Fig. |U we show the on-top correlation hole for the ex- 
act RPA, the RPA-based LSDA (see Ref. [H) and the 
RPA-based GEA of Eq. JT2JI. We see that for a jel- 
lium surface the RPA-based LSDA on-top correlation 
hole nearly coincides with the corresponding exact RPA 
on-top correlation hole; this is in contrast with the case 
of strong inhomogeneous systems (e.g., Hooke's atom)j^ 
The gradient correction of Eq. (fl"2"|) improves the already 
accurate RPA-based LSDA on-top correlation hole in the 
slowly-varying density region, but is inacurate in the tail 
of the density. Fig. 2] also shows that the integrated 
b c (z, z) of Eq. (fTTj) is more (less) negative in the vacuum 
(bulk) than the actual on-top correlation hole. 

At this point, we would like to emphasize that while 
the RPA on-top correlation hole in the bulk is too nega- 
tive but finite, the on-top correlation hole diverges in the 
bulk within a TDDFT scheme that uses a wavevector 
and frequency independet xc kernel like in the adiabatic 
local-density approximation (ALDA) 



f£? A [n](T,x>, U ) = 



dn(r) 



(13) 



or the energy-optimized local-density approximation of 
Ref. H^. (See the discussion after Eq. (3.9) of Ref. \x 



Here, v^ nl} [n(r)] is the xc potential of a uniform 
electron gas of density n(r). An xc kernel borrowed 
from a uniform-gas xc kernel that has the correct large- 
wavevector behavior (see, e.g., the xc kernels of Refs. |4G. 
I4HI42I ) would yield a finite on-top correlation hole. Fig. 



[5] shows the integrated correlation hole of Eq. (fPTj) for an 
electron at the vacuum side of the surface, at the posi- 
tion z = — 0.5\p and for several values of the electron- 
density parameter r s : 1.5, 2.07, 3, 4, 5, and 6. In the 
bulk, the correlation hole exhibits damped oscillations 
with r s -dependent amplitude and a period that does not 
depend on the electron density and is close to the period 
(~ 0.56Aj?) of the corresponding oscillations exhibited by 
the exchange-only hole. 

Finally, we look at the xc energy density e xc de- 
fined in Eq. ([3]). We note that adding to the ac- 
tual e xc of Eq. ^ an arbitrary function of the posi- 
tion r that integrates to zero yields the same total xc 
energy The Laplacian of the density V 2 n integrates to 
zero for finite systems, it plays an important role in the 
gradient expansion of the kinetic- energy densit y 44 ' 45 ! 46 
and it is an important ingredient in the construction of 
density-functional approximations for the kinetic energy 
densit y 44 ! 45 and the xc energy.— 

We define the simplest possible Laplacian-level RPA- 
based LSDA (the RPA-based L-LSDA) xc energy density: 



^L — LSDA— RPA 



LSDA— RPA 



(r) -CV 2 n(r), (14) 



where C is a constant parameter which we find by mini- 
mizing the difference between e x ^ A ~ L ~ LSDA and e^f . 
We find C = 0.3 for a jellium slab with r s — 2.07, and 
its value gets larger as r s increases. 

In Fig. M we show Ae xc (z) = e^f A {z) - e a vp rox (z) 
versus z/Xp for a jellium slab with r s = 2.07 and several 
RPA-based approximations for e a x p c prox {z) . The RPA- 
based PBEi 5 " improves considerably the behavior of the 
RPA-based LDA. The ARPA-GGAil is a GGA func- 
tional that fits the RPA xc energy density of the Airy 
gas and is remarkably accurate for jellium surfaces. The 
RPA-based GGA++ is the RPA version of the GGA++ 



of Ref. 
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FIG. 6: Ae xc (z) = e« PA (z) - e° 



versus z/Xf at a 



surface of a jellium slab, for several xc approximations: RPA- 
based LSDAA RPA-based PBFA ARPA GGA^, RPA- 
based GGA+4A and RPA-based L-LSDA (Eq. $T4$ with 
C — 0.3). The bulk parameter is r s — 2.07, and the edge of 
the positive background is at z = 0. 



fined in Eq. (3) of Ref. \M-) Although the GGA++ func- 
tional was constructed for the Si crystal, we observe that 
the RPA-based GGA++ improves over the RPA-based 
LSDA in the bulk near the jellium surface showing that 
it can be a good approximation for systems with small os- 
cillations. (In the bulk, close to the jellium surface, there 
are Friedel oscillations as well as quantum oscillations 
due to the finite thickness of the jellium slab). We note 
finally that e^ A ~ L ~ LSDA significantly reduces the local 
error of the RPA-based LSDA near the jellium surface, al- 
though by construction E*f A - L - LSDA = E PPA ~ LSDA . 



density at a metal surface. When the electron is in the 
vacuum, its hole remains localized near the surface (its 
minimum is on the image plane) and has damped oscilla- 
tions in the bulk. We find that the on-top correlation hole 
is accurately described by local and semilocal density- 
functional approximations, as expected from Ref. HI . We 
also find that for an electron that is localized far outside 
the surface the main part of the corresponding xc hole 
is completely localized at a plane parallel to the surface, 
which is the image plane. 

Because of an integration by parts that occurs in the 
underlying gradient expansion, a GGA (or meta-GGA) 
hole is meaningful only after averaging over the elec- 
tron density n(r) i 18 i 20 This average smooths the sharp 
cutoffs used in the construction of the angle-averaged 
GGA xc hole density. The wavevector analysis of the 
jelium xc surface energy is an important and hard test for 
the LSDA, GGA, and meta-GGA angle-averaged xc hole 
densities, showing not only the accuracy of the xc hole 
but also the error cancellation between their exchange 
and correlation contributions. Thus, Refs. HH an d 20 i 21 
have shown that the TPSS meta-GGA 2 ! and the PBEsol 
GGA3S xc hole densities improve considerably the ac- 
curacy of their LSDA and PBE counterparts at jellium 
surfaces, both within RPA and beyond RPA4& 

The exchange energy density does not have a gradient 
expansion^, as does the kinetic energy density. However 
the existence of gradient expansion of the xc energy den- 
sity is still an open problem. We use our RPA xc hole 
density to compare the xc energy densities of several ap- 
proximations. The most accurate ones are ARPA GGA 
of Ref. E3 and RPA-based L-LSDA of Eq. CO 
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III. CONCLUSIONS 

We have presented extensive self-consistent calcula- 
tions of the exarf-exchange hole and the RPA xc hole 
at a jellium surface. 

We have presented a detailed study of the RPA xc hole 
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